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We investigate the thermoelectric properties of a system near a pressure driven Mott-Hubbard 
transition. The dependence of the thermopower and the figure of merit on pressure and temperature 
within a degenerate Hubbard model for integer filling n = 1 is calculated using dynamical mean 
field theory (DMFT). Quantum Monte Carlo method is used to solve the impurity model. Obtained 
results can qualitatively explain thermoelectric properties of various strongly correlated materials. 
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The discovery new strongly correlated materials and 
improvement of theoretical methods offer new perspec- 
tives in the search for systems with good thermoelec- 
tric performance jjj. The dimensionless figure of merit, 
ZT = S 2 aT ' I k, provides us with a quantified measure of 
the thermoelectric performance of materials. The higher 
values of the thermoelectric figure of merit corresponds 
to the better thermoelectric properties. To get maximal 
thermoelectric response at a fixed temperature, T, we 
need to have maximum possible thermopower, S, and the 
electrical conductivity, cr, while the thermal conductivity, 
k, should be the smallest possible. The thermal conduc- 
tivity has two contributions electronic, n e , and lattice, 
kl, one, n = n e + kl. For some times, the highest value 
of the figure of merit was equal to one and only a little 
improvement in getting materials with higher figure of 
merit has been achieved in the last two decades. Since 
the electronic structure of strongly correlated electron 
systems exhibit properties which have no analogy with 
those of weakly correlated compounds, it is important to 
understand their physical properties and how they im- 
pact the thermoelectric power. Strongly correlated elec- 
trons are sensitive to small changes in their control pa- 
rameters such as temperature, doping and pressure and 
hence require a detailed investigation. There is now a 
strong interest in understanding how this changes affect 
the thermoelectric properties in order to find more effi- 
cient material with higher thermoelectric response. Even 
a modest increase in ZT could substantially impact a 
number of applications jlj . 

Two important parameters in strongly correlated elec- 
tron systems are the carrier concentration and the ratio 
of the on-site interaction U to the bandwidth W. This 
ratio can be altered by applying external pressure or in- 
ternal pressure by means of isovalent substitutions. 

In this paper we study behaviour of the figure of 
merit and the thermopower near a pressure driven Metal- 
Insulator Transition (MIT) within framework of the de- 
generate Hubbard model. Our focus is on the effects 
of orbital degeneracy. There are several motivations for 
our study. First, in many strongly correlated materials 
the orbitals, which bring a major contribution into un- 
usual physical properties, are degenerate. Second, the 



orbital degeneracy allows us to study the effects of par- 
ticle hole asymmetry in the Mott insulating state with 
filling n — I in a natural way (one can also introduce 
particle hole asymmetry in this case by adding hopping 
integrals beyond nearest neighbors). Third, the degen- 
eracy will allow us to use the integer filling n = I what 
in the real system corresponds the situation when the 
number of the electrons is equal to the number of sites. 
We use Coulomb repulsion U and temperature T as the 
parameters to vary and take the half bandwidth D to be 
unity. 

While the calculation of the numerical value of the 
thermoelectric coefficient requires the detailed modeling 
of the band structure and interaction constants of the 
relevant compounds, we expect that the qualitative fea- 
tures related to the thermoelectricity of a Mott insulator 
at integer filling can be captured by the simplest model, 
the degenerate Hubbard model, and we undertake this 
study in the paper. Indeed the trends discussed here, 
are found in NiS2- x Se x and Nii_ x Co x S2 as we discuss 
toward the end of this paper. The algorithms necessary 
for more realistic calculations of thermoelectricity will be 
discussed in a future publication Q . 

To treat the two-band degenerate Hubbard model we 
use the Dynamical Mean Field Theory (DMFT) §|. Re- 
cent development of the DMFT has given a boost to the 
study of strongly correlated systems and their properties. 
In particular it has resulted in a detailed understanding 
of the Mott transition. 

Previous studies of thermoelectric properties using 
DMFT were done within framework of the one-band 
Hubbard model strong temperature and doping 

dependence of the thermopower was reported. Palsson 
and Kotliar considered the orbital degeneracy away from 
half filling in the limit of infinite interaction strength . 
A good description of the Seebeck coefficient dependence 
on temperature and doping for titanates was obtained. 

In this paper we focus on the integer occupancy case 
with orbital degeneracy. We also keep U finite and study 
dependence of thermoelectric properties on the interac- 
tion. Changing the ratio of the interaction strength one 
can mimic the effects of pressure. 

The TV-fold degenerate Hubbard Hamiltonian reads: 
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H = - iyct^Cjv + — ^ n ia n i<T > - p^n^, (1) 

i.a^a 1 La 

where (ij) runs over nearest neighbor sites, a is the spin 
and orbital index which runs from 1 to N . The two- 
band degenerate Hubbard model corresponds to N = 4. 
The hopping matrix is given by ty, U is the Coulomb 
repulsion and p is the chemical potential. We also want 
to restrict ourselves to the paramagnetic phase (in both 
spin and orbit indices) i.e. we are interesting in metal- 
insulator transition and its influence on transport prop- 
erties due to electronic correlations only, ft is worth to 
mention that orbital degeneracy suppresses the magnetic 
correlations enough to make our assumption easily ful- 
filled j§] . To make computations easier further we assume 
t^j t. 

In this paper we use the model on a cubic lat- 
tice. For calculation of thermodynamic properties we 
use the semicircular density of states (DOS) p(e) = 
(2/tt£») v /1 - {e/D) 2 , with the half-bandwidth D = 
W/2 — 2t. The semicircular DOS corresponds to an in- 
finite coordination Bcthc lattice which, as it was found 
earlier |3) , gives a good description of three dimensional 
systems. 

The next step is to solve the Hubbard model. The stan- 
dard way in DMFT to do it is to map the lattice Hub- 
bard model onto the effective impurity problem, which 
is a generalized single impurity Anderson model, where 
the operators carry an orbital index. Supplemented with 
the self-consistency condition jjj solution of the impu- 
rity model gives us the solution of the original Hubbard 
model. Quantum Monte Carlo method with extended 
Hirsch-Fye algorithm [9[JlO|] is used as the impurity solver. 
To reduce computational errors (the Trotter breakup) 
imaginary time interval is set to At — 1/4. During 
the self-consistency procedure one needs to make direct 
and inverse Fourier transformations for the Green's func- 
tions (GF). In the present calculations we use a modified 
Fourier transformation (FT) in order to get correct re- 
sults (for details see Appendix). 

The output of the self-consistent procedure described 
above is the Green's function on imaginary time axis. To 
calculate the thermoelectric properties we need to know 
behaviour of GF on real frequency axis. The Maximum 
Entropy (ME) method is used to make the analytical 
continuation of imaginary-time GF to DOS on the real 
frequency axis. Knowing imaginary part of GF we re- 
construct frequency dependence of the real part G(u>) 
through Kramers-Kronig relations. Obtained GF is used 
to calculate the transport coefficients. 

Using Kubo formalism one can express thermoelectric 
coefficients in terms of current-current correlation func- 
tions which are reduced within DMFT to averages over 
the spectral density function p(e,w) J7|,pd|: 
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where the coefficients A n have the following form, 

where T is temperature, /(w) is the Fermi distribution 
function, s describes spin and orbital indexes which runs 
from 1 to Ndeg- The spectral density, p s (k,u>), and the 
transport function, $(e), contain the relevant informa- 
tion about the bare band structure, e^: 

•<«>- (4) 

As it is seen from the above formulae (||) and (|4|) contri- 
bution to the thermopower strongly depends on tempera- 
ture. For low temperatures only states close to Fermi sur- 
face (FS) contribute to thermoelectric properties while 
for high temperatures the entire Brillouin zone is im- 
portant. Coulomb repulsion, U, acts on the transport 
coefficients via changes in the spectral density function, 
p(e,cu). There is a very simple mnemonic rule to define 
the sign of the thermopower. For small temperatures it 
depends on DOS slope at the Fermi energy: if the DOS 
curve rises when it crosses the Fermi energy then the sign 
is negative and vice versa. For very high temperatures 
it depends on weight of the Hubbard bands: if weight 
of the Hubbard band above the Fermi energy is larger 
than contribution from the band below the Fermi level 
then the thermopower is negative. This rule will help us 
easily understand the thermopower behaviour analyzing 
temperature dependence of DOS. 

In Fig. ^ we plot dependence of the thermopower on 
interaction strength U. The thermopower changes sign 
in a region of Coulomb repulsion 2.5 < U < 3, the region 
where the system undergoes metal-insulator transition at 
low temperatures. The MIT for high temperatures trans- 
forms into the metal-insulator crossover what is reflected 
in the spectral function (DOS) behaviour presented in 
Fig. 0. It is worth to notice quite strong S dependence 
on U in that crossover region. When U — > the ther- 
mopower will have a finite value depending on tempera- 
ture and in opposite limit U — > oo we expect saturation 
of the thermopower dependence which will be also tem- 
perature dependent. The interaction dependence of the 
thermopower helps us qualitatively understand the pres- 
sure influence onto the system as well as behaviour of 
a system (e.g. NiS2_ x Se x ) changing its properties from 
insulating to metallic ones by varying the bandwidth. It 
follows from the fact that the only important parameter 
in the system is the the ratio U/W and one can vary any 
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of two variables to get qualitative description of a system 
properties. 

Behaviour of the figure of merit on interaction strength 
(Fig. ||) is possible to understand from its definition. It 
has minimum where the thermopower changes sign. Af- 
ter the sign change it substantially increases and has 
maximum for higher values of U. Rather high values 
of ZT (more than one) can be explained by the absence 
of the lattice contribution in the thermal conductivity. If 
one takes into account the lattice thermal conductivity 
the figure of merit drops to value below one. Usual val- 
ues of lattice contribution into the thermal conductivity 
in transition metal-oxides lay in region 0.1-10 W/mK. We 
believe that 1 W/mK is a reasonable value for pyrites. 
We see that the figure of merit is maximal for interaction 
U between 4 and 4.5, but even in this region obtained 
figure of merit could not be competitive with the one for 
semiconductors due to the lattice contribution into the 
thermal conductivity. 

Temperature dependence of the thermopower for U = 
3 presented on Fig. ^ can be easily explained from anal- 
ysis of DOS temperature dependence in Fig. ^[ For very 
low temperatures the thermopower has linear dependence 
(one can show it analytically, it is nearly impossible at 
the present stage to reach very low temperatures with 
QMC) and negative sign indicating presence of quasi- 
particles in the system what is clearly seen in Fig. |B| 
(/3 = 16) where a good quasiparticle peak at the Fermi 
energy is observable. With increasing temperature a 
smooth crossover from quasiparticles excitations to in- 
coherent ones reflected in non-monotonic thermopower 
dependence which changes sign from negative to posi- 
tive one. In Fig. ^| for temperatures /3 = 12 and (3 — 8 
we see how quasiparticle peak of DOS gets lower and 
then disappears. From this observation we can conclude 
that positive value of S tells us that the band structure 
consists from two Hubbard subbands only. Similar be- 
haviour of the thermopower (with the sign change) occurs 
in the single band Hubbard model upon doping ||ll[j and 
in the Hubbard model with frustration p] as well as in the 
periodic Anderson model |jl2| . Further temperature be- 
haviour of the thermopower is new one. At temperature 
(3 = 8 in Fig. || the thermopower reaches its maximum 
and then starts decreasing (we should notice here that 
temperature /3 — 4 is the coherent temperature for the 
studied system) . For (3 = 4 the spectral function consists 
from two bands, where lower Hubbard band, which is 
closer to the Fermi energy has greater contribution into 
the thermopower than the upper one, which is located a 
bit further from the chemical potential position. With in- 
creasing temperature two bands became less asymmetric 
to the Fermi energy and the weight from the upper Hub- 
bard band (n = 1) is larger than the contribution from 
the lower band. It means that the thermopower should 
became negative as it seen in Fig. [| (the last point). And 
finally, for very high temperature the two bands collapse 



into the one. But the sign of the thermopower remains 
unchanged by the same reason as we discussed above. 
The thermopower sign change in the high temperature 
limit due to Mott-Hubbard bands collapse was reported 
in 0. 

To gain further understanding of the temperature and 
interaction dependence of the thermopower we study its 
high- frequency behaviour in the high-temperature limit. 
To do so we generalize the thermoelectric response to 
finite frequencies (ac-thermopower) [pTj: 
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where coefficients L 11 and L 12 are defined as 
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Here e~^ a = Tr (e - ^-"^), i,j = x, J x x = j x , 
J 2 = jx I Jx and j® are electrical and heat currents in x 
direction correspondingly, q x is the polarization operator 
satisfying j x (t) = dq x (t)/dt. 

Expanding numerator and denominator of the ther- 
mopower (L n (w) and L 12 (ui)) in frequency and dividing 
one onto another we obtain the following expansion of 
the thermopower in the high-frequency limit: 
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The relevant commutators are given by: 
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This expression consists of three terms. The first 
one, proportional to t, doesn't contribute to the high- 
temperature expansion if one takes into account only 
nearest neighbour hoppings. The second term (propor- 
tional to chemical potential) similar to what we have 
in the denominator and the last one proportional to U 
contribute to high-temperature expansion of the ther- 
mopower. The final result we end up with is presented 
below 
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where n corresponds to contribution to filling n coming 
from Ndeg — 1 degrees of freedom (spin and orbital) and 
is equal to a i^ a ((l — n a )n<T'}- If one assumes high- 
temperature behaviour of the chemical potential in the 



3 



form fi/T — a, where a should be negative for the filling 
n < T^Ndeg, one obtains that the sign of the thermopower 
should be negative in the high-temperature limit. Hence 
ac- and dc- thermopowers have the same sign in the high- 
temperature limit. 

It is quite difficult to make a direct comparison of a 
model calculation and situation in real materials due to 
complexity and wide variety of real structures. To make a 
comparison with our model calculations we need a mate- 
rial with degenerate e g band hosting one electron (hole). 
As we mentioned in the introduction a good candidates 
for the comparison are pyrites compounds with doubly 
degenerate 3d — e g Ni band. An attempt to understand 
high-temperature dependence of the thermopower in this 
materials on base of one-band electron correlation theory 
with a good fit of experimental data using six parame- 
ters model was given some years ago in work of Kwizera 
et al. |15|. To explain low temperature dependence of 
the thermopower it was suggested that both holes and 
electrons participate in charge transport and a two-band 
model was proposed as an appropriate model for data 
interpretation jp|] . Our calculations show that the two- 
band degenerate Hubbard model treated within DMFT, 
is sufficient model for at least qualitative description of 
pyrites in the whole temperature range [0;[^,[16| . In the 



case of half-filled Ni 3d — e g band in NiS2- x Se x theoret- 
ical model under consideration due to the symmetry in 
the system should give zero thermopower for all tempera- 
tures, while experimentally it is zero only for temperature 
T < 100 K. For higher temperatures the thermopower 
becomes positive. It is clear that one needs to take into 
account effects of other bands (filled ti g band in the first 
turn). The thermopower behaviour in Nii_ x Co x S2 with 
temperature looks rather similar as in NiS2- x Se x but 
with Co substitution of Ni atoms occupancy of 3c? — e g 
band changes from two (x = 0) to one {x = 1). One can 
track the experimental situation in NiS2- x Se x analyzing 
model calculations (for temperatures low enough to ne- 
glect contributions coming from other bands) especially 
in the case C0S2 which corresponds to one-quarter filled 
3d — e g band. 

In conclusion, we calculated dependencies of the ther- 
mopower and the figure of merit on interaction strength 
and temperature of the thermopower in the two-band de- 
generate Hubbard model for integer filling n = 1. The 
strong dependence on studied parameters was obtained. 
Analytical estimations of the high-frequency limit of ac- 
thermopower is provided. In the high-temperature limit 
ac- and dc- thermopowers have the same sign (negative 
for filling n < 0.5). We expect that the thermopower 
behaviour in the three-band degenerate Hubbard model 
would behave in a similar way. It means that we can 
qualitatively understand behaviour of the thermoelectric 
properties in a wide range of strongly correlated materi- 
als where a major role played by d- and /-shell electrons. 

Systems near the temperature driven Mott transi- 



tion exhibit very rich thermoelectric behaviour. As we 
cross the localization-delocalization threshold the ther- 
mopower increases substantially. There are the high- 
temperature precursors of the first order Mott transi- 
tion that takes place at lower temperature. As the tem- 
perature is lowered, and the critical Mott endpoint is 
approached, the specific heat diverges and the entropy 
jumps. The figure of merit is low in this crossover regime 
because of the great deal of cancellation between the 
quasiparticle contributions which are in our model elec- 
tron like, and the Hubbard band contributions which are 
hole like. A detailed phase diagram temperature versus 
filling for different values of U using the current approach 
will be studied elsewhere. 
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APPENDIX A: FOURIER TRANSFORMATION 
IN QMC CALCULATIONS 

In the self-consistent procedure to solve the impurity 
problem using QMC method j| we need to do two di- 
rect Fourier transformations for GFs G(t) and Gq(t) 
(the Weiss function) and one inverse Fourier transfor- 
mation (iFT) for Gq(t) only. Function Gq(t) is an in- 
put for QMC simulations producing G(t). Two direct 
Fourier transformations are necessary to impose the self- 
consistency condition which is usually written in fre- 
quency space. Both, GF and the Weiss function contain 
the discontinuities at r = and t = (3. The disconti- 
nuities in these functions and their derivatives determine 
the high-frequency behaviour of their frequency depen- 
dent Fourier transformations G(ui) and Go(ui). As it is 
well known that information about high-frequency be- 
haviour of GFs is absent in QMC simulations themselves 
(maximum frequency available is the Nyquist frequency 
^max = 1/2 A, where is A is imaginary time interval). 
The high frequency (small imaginary time) information 
is not contained in the QMC itself but has to be incor- 
porated into the splining procedure using additional in- 
formation. 

This information is available from calculation of corre- 
sponding moments of the GFs. To make this connec- 
tion clear, we make consecutive integration by parts of 
the Fourier integral: 



f (-1)W(GW(0)+GW^)) , vn 
G{iuJ n ) = pr ; ^ fc+1 (Al) 

fc=0 



(iui n ) k 



(-1) 



N+l 



N+l 



.d N+1 G{r) 



dr. 



Values of GF sum and it's derivatives can be expressed 
via corresponding momenta. These momenta can be cal- 
culated similar to the one-band Hubbard model jf?]]. A 
momentum of k degree is defined as follows: 



M (fe) = / du>ui k p(ui). 



(A2) 



We can bind Eqs. (Al) and (A2) writing the follow- 
ing expression for sum of GFs and it's derivatives in 
imaginary-time space: 



G( fc) (0) + G (fc) (/3) =Af (fc) , 



(A3) 



where k — 0, . . . N . To make direct Fourier transforma- 
tion, first, we interpolate G(r) defined in L points (L- 
number of time slices) and then take analytical Fourier 
transformation of resulting function. The procedure of 
interpolation plays a key role. We used the cubic spline 
interpolation for G(r) where condition of continuous sec- 
ond derivative (G^) is imposed. Using this condition we 
write a system of linear equations to find interpolation 
coefficients. To close the set of linear equations we need 
to put correct boundary conditions. The standard way to 
do it is to set the second derivatives at the end points of 
the time interval [0,(3] to zero (so called, natural spline). 
In our approach we use analytical information about the 
momenta (sum of the second and first derivatives at the 
end points, M (1) and M^). 

To obtain this information we need to know expansion 
of the self-energy in 1/w series or in another words we 
need to know first two momenta of the self-energy. 



il4 0) = {2N - \)Un, 
r(i) 



(A4) 



M£> = (2N - l)U 2 n(l - n) + 2U 2 C% N - 1 (nn), 

(fe)l 



where N is the number of bands and Cz — U Y'' M is a 

n n\{k—n)\ 

combinatorial factor which arrives due to the spin and 
orbital degeneracy. 

The self-energy expansion contains density-density cor- 
relations functions for different spins and orbitals. One 
way to obtain them is to use an approximate scheme 
which is accurate at high frequencies such as the coher- 
ent potential approximation (CPA) jnl. Another 
possibility is to evaluate the correlation functions which 
enter in Eq. (A4) using the QMC procedure which allows 
us to compute arbitrary local correlation functions such 
as the density-density correlators. The self-consistency 
procedure ensures us that they are correct correlation 
functions when the self-consistency is reached. Having 
correct momenta we solve the set of linear equations to 
obtain the cubic spline function which is Fourier trans- 
formed analytically. 

The power of modern computers allows us to use in 
imaginary time (r) space up to L — 256 time slices only, 
while in frequency (uj) space we do not have such limita- 
tions and can use as much frequency points as we want 
but only frequencies less than the Nyquist frequency have 
physical meaning (usually it is less than 2 8 ). As we know 
GF has 1/uj asymptotic behaviour and dealing with the 
finite number of frequency points we introduce the finite 
energy cut-off. Cutting off the tail of GF, we remove the 
discontinuity from the Fourier transformed function in 
t space. To correct the situation we subtract the high- 
frequency tail from GF and make FT of the obtained 
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function numerically and make Fourier transform of the 
tail analytically. Finally, in t space we sum up the ob- 
tained functions and as the result we have the correct 
inverse Fourier transformation, we call it the "new" one. 

To demonstrate the difference between two "old" and 
"new" Fourier transformations we calculate the self- 
energy in the considered model for U = 4 and (3 = 8. 
Results are plotted in Fig. || where the "new" solution 
is drawn by solid line and the "old" one is plotted by 
dashed line. We see that the "old" one has the region 
where the self-energy changes sign ("overshoots"). It cor- 
responds to unphysical contribution to self-energy which 
should always keep the same sign for positive or negative 
frequencies. 

Finally, we stress that difference between two Fourier 
transformations becomes substantial especially in critical 
regions of parameters: low doping, high values of U and 
low temperatures T. 



FIG. 1. Dependence of the thermopower, S ( in units of 
fcs/e = 86/j.V/K ), on Coulomb interaction, U (in D = 1 eV) 
for (3 = 8. 



FIG. 2. Spectral functions (DOS) for U = 2.5 and 3 for 

/3 = 8. 

FIG. 3. The figure of merit, ZT, vs. interaction strength, 
U, for temperature /3 = 8 and different values of the lattice 
conductivity kl [W/mK]. 

FIG. 4. Temperature dependence of thermopower, S for 
U = 3. The solid line is guide for the eye. 

FIG. 5. Temperature dependence of density of states for 
U = 3. 

FIG. 6. Dependence of imaginary part of the self-energy on 
Matsubara frequency axis to (/3 = 8) for two different Fourier 
transformations 'old' (dashed line) and 'new' (solid line) one. 
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